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Abstract 

The effect of velocity correlations on the equal-time density autocorrelation function, e.g. the 
pair distribution function or pdf, of a hard-sphere fluid undergoing shear flow is investigated. 
The pdf at contact is calculated within the Enskog approximation and is shown to be in good 
agreement with molecular dynamics simulations for shear rates below the shear-induced ordering 
transition. These calculations are used to construct a nonequilibrium generalised mean spherical 
approximation for the pdf at finite separations which is also found to agree well with the simulation 
data. 
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I. INTRODUCTION 



In equilibrium, simple fluids exhibit spatial correlations which are characterized by the 
pair distribution function (pdf) describing the probability of finding two atoms with a give 
relative orientation and separation. Equilibrium liquid state theory is primarily concerned 
with the calculation of the pdf and a number of successful approaches have been developed 
including the Percus-Yevik approximation for hard-spheres, the mean-spherical approxi- 
mation and the more recent self-consistent integral equations]]]]. Knowledge of the pdf is 
equivalent to knowledge of the density-density static correlation function |IJ and once this is 
known, all other interesting static correlation functions, e.g. density-energy, energy-energy, 
are immediately known because the velocity-dependence of the two-body distribution 
function in equilibrium is trivial. It is a characteristic of nonequilibrium fluids that this 
property no longer holds [0, |3|, §,0], and the presence of velocity correlations is the reason 
that the determination of static correlations in nonequilibrium fluids, over all densities and 
length scales, is a difficult problem. 

The velocity correlations that occur in nonequilibrium fluids, as well as in fluctuations 
about the equilibrium state, are generated by collisions which have the effect of altering 
the two-body probability distribution so that even if the velocities of the atoms prior to 
a collision are assumed to be independent variables, the velocities after a collision are not 
independent. The question of whether the of velocities of two atoms prior to a collision are 
really independent variables has been much studied in statistical mechanics over the last 30 
years and phenomena such as long-time tails and long-ranged correlations are proof that this 
assumption is not strictly adhered to|| [7|, although in many cases it remains a good approx- 
imation. While the calculation of the pre-collisional correlations is a very difficult problem, 
it has recently been noted that, for fluids interacting via a hard-core potential, it is possi- 
ble to calculate the post-collisional correlations in an arbitrary nonequilibrium state for the 
special case of the two atoms being in contact using the same assumptions as underlie the 
Enskog theory of the one-body distribution function. This allows one to calculate all static 
correlation functions for two atoms in contact up to this level of approximation. It was 
subsequently shown that this information could be combined with a formalism borrowed 
from equilibrium liquid-state theory to create a successful model of the pair distribution 
function of a granular fluidQ (i.e., a fluid of inelastic hard-spheres). The purpose of the 
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present paper is to describe an extension of this model to inhomogeneous systems and to 
examine its application to the particular case of a fluid of elastic hard-spheres undergoing 
uniform shear flow (USF) and to present detailed comparisons of the theory to the result 
of molecular dynamics simulations. Uniform shear flow, in which the velocity in one Carte- 
sian direction varies linearly with position along another axis, is a particularly interesting 
example since the density-density correlation function can be studied experimentally by 
means of light scattering ||. Furthermore, the hard-sphere model is generally accepted as 
a reasonable analogy to certain types of colloidal suspensions, see for example ref. || and 
references therein, for which it is possible to achieve conditions of strong shear (e.g., shear 
rates comparable to the mean free time of the colloidal particles) in the laboratory which 
are otherwise inaccessible in simple fluids. 

The second section of this paper reviews the theory behind the calculation of static cor- 
relations at contact and evaluates the density-density correlation function at contact for the 
special case of USF. This makes use of recent work on solution of the Enskog equation for 



high shear rates [|T0|, PH to extend an earlier calculation resulting in an explicit expression 
for the pdf at contact in a sheared fluid. The third section deals with models for the pdf 
at finite separations. It reviews two well known theories which are potentially applicable to 
atomic length-scales: the kinetic model studied by Hess and co-workers [TI| and the Langevin 
model of Ronis[|T^]. The former involves an undetermined parameter which, if the theory 
is to apply to atomic length scales, can now be fixed by requiring agreement with the cal- 
culations for two atoms in contact. The latter, while not involving any a priori unknown 
parameters is nevertheless phenomenological and a diffusion constant appearing in its for- 
mulation has in fact been treated as a free parameter when comparing to experiment [§]. 
Again, it is noted that the parameter can be fixed unambiguously by requiring agreement 
with the calculated value at contact. It is also shown that these two theories are in fact 
very closely related not withstanding their different motivations. Finally, in this section 
the nonequilibrium version of the Generalized Mean Spherical Approximation (GMSA) is 
introduced as a means of modeling the pdf at finite separations based on the atomic-length 
scale information coming from the calculations of Section 2 and, qualitatively, the large- 
separation (i.e., small wave- vector) information provided by mode-coupling theories, based 
on either kinetic theory |T2], |TJ], [T5|] or fluctuating hydrodynamicfTB|, |TE| , |T7| , 113, of which the 
Ronis theory is an example. This is not unlike the original motivation of Weisman in intro- 
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during the equilibrium GMSA as a means of improving on the Percus-Yevik approximation 
by incorporating accurate knowledge about the pdf at contact, from the Carnahan-Starling 
equation of state and the pressure equation, to construct a model of the pdf for an equilib- 
rium hard-core fluid accurate over a wide range of densities [IS]. Recent work by Yustes and 
Santos [EO, I2TI, E2fl, as well as Carraro and Ciccariello[E3|, has shown first that the Percus- 



Yevik approximation may be viewed as, in some sense, the simplest approximation that 
provides certain analytic properties that any distribution function must satisfy and second, 
that the GMSA of Weisman may be viewed as a framework for systematically extending 
this model so as to incorporate additional constraints. It is with this motivation that the 
extension of the GMSA to nonequilibrium systems was proposed^ as a means of modeling 
the density-density correlation function, i.e. the pdf, at all length scales. 

In section 4, these calculations are compared with the results of molecular dynamics 
simulations over a wide range of shear rates and densities. As noted previously]^], there 
seems to be a strong correlation between the rapid decrease, with increasing shear rate, of 
the pdf in certain directions and the onset of shear-induced ordering of the fluid. Below 
this transition, it is shown that the Enskog calculations of the pdf at contact are quite 
accurate at small shear rates and low densities and becomes increasingly inaccurate as the 
density and/or the shear rate increases. The theories for the pdf at finite separations are 
also compared to MD and it is found that all three theories are in qualitative agreement 
with the GMSA providing the best quantitative agreement with simulation. The paper ends 
with a discussion of the prospects to extend these results to other systems. A preliminary 
description of some of these results has appeared previously ||. 



II. THEORY OF CORRELATIONS AT CONTACT 



A. Hard-sphere statistical mechanics 

Consider a system of N elastic hard spheres of diameter a in a cubic volume V = L 3 
described by a Cartesian coordinate system with axes x, y and z. The boundary conditions 
will be discussed below. The dynamics of the atoms consist of free-streaming, subject to the 
boundary conditions, interrupted by elastic collisions. Two atoms having coordinates q,, p« 
and qj,Pj at time t will collide at time £_ provided that a = |qy(t_)|, where qr,(t_) = 
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qi(t-) — <ij{t~) and provided q^(i_) • Py(i_) < 0. Immediately after the elastic collision, 
the momenta become 



Pi(t+) = p;(t-) - % (t_) (%(*-) • Py(*-)) 
Pj(*+) = Pi(*-) + ■ Py(*-)) 



(1) 



so that the relative momentum is reversed along the line of contact and the total momentum 
is unaffected. 

The statistical description of the system is characterized by the N-body distribution, 
Pn(%i, X2---XN', t) which gives the probability of finding the system at a given phase point, 
where atom 1 has phase X\ = (qi,Pi) etc., at the specified time t. Its evolution is specified 
by the pseudo-Liouville equation 

N „ N 



where the collision operator is given byp4 



p N {x 1 ,x 2 ..-x N ;t) = 



(2) 



T-(ij) 



dqij S {qij - a) 



bij - 1 



(3) 



with the effect of the momentum transfer operator by on an an arbitrary function is to replace 
the relative momentum p^- by its post-collisional value — 2qjj (q^- ■ py) (see Eq. |]). The 
final term of Eq. (0) describes any external one-body forces acting on the atoms. Integrating 
(§) over N — n of the coordinates yields the n-the equation of the BBGKY hierarchy which 
the relates the n-body distribution to the n + 1-body distribution. In particular the result 
of choosing n = N — 1 is 



d 



dt <9qi ' dpi 



Pi{xi,t) 



(4) 



(N - 1) / dq 12 S (q 12 - a) 



'12 



@ (~qi2 • Pl2)P2(^l,^2;0" 



If the two-body distribution on the right is approximated by pzixu %2', t) ~ 
Pi(xi,t)pi(x2',t)g(qx,q2',t), the result is the well-known Enskog equation for the one-body 
distribution of a system of hard spheres (here, p(qi,q2;t) is the probability to find two 
atoms at positions qi and q 2 and is normally approximated by the equivalent local equi- 
librium function). In fact, examination of eq.(|4]) shows that the necessary approximation is 
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actually 

5 (gi 2 - cr) (-qi2 • P12) P2(xi,x 2 ] t)xi6 (qu - a) (-qi 2 • P12) Opi( x 2; *)#o(qi, q2; t) 

(5) 

which is a somewhat weaker approximation than the assumption that the two-body distri- 
bution always factorizes: rather, one need only assume that it factorizes for the case that 
the two atoms are in contact and approaching one another which is to say, just prior to 
a collision. This is a precise statement, for hard-spheres, of Boltzmann's "assumption of 
molecular chaos" . Immediately after a collision, the direction of the relative momentum is 
reversed and the momenta of the two atoms are obviously correlated. In fact, it has been 
shown[Q that the approximation given in Eq.(^) implies the form of the entire two-body 
distribution at contact is given by 

5 (qu - cr) p 2 (x 1 ,x 2 ;t) ~ 6 (q 12 - cr) p 1 (x 1 ;t)p 1 (x2;t)g (q 1 ,q 2 ;t) 

+$ 0?i2 - <?) © (qia • P12) 612-I Pi{xi]t)pi{x 2 ]t)g {qi,q 2l t{e) 

The distribution is seen to have to parts: the first term on the right describing uncorrelated 
atoms the second term describing velocity correlations which arise because of collisions. In 
equilibrium, the second term vanishes but for non-equilibrium systems, it is generally present 
and can give rise to substantial structural effects as will be discussed below. This relation is 
critical in that it can be used to calculate, to the same level of approximation as is inherent 
in the Enskog equation, any static two-body correlation function at contact. 

Finally, the nonequilibrium pair distribution function is defined, as in equilibrium, by 

0(01. Qa;*) = v J dp 1 dp2p 2 (x 1 ,x 2 ',t). (7) 

From the definition of the local density field 

N 

n(r) = ^5(r-qj (8) 

it follows that the pdf is related to the density autocorrelation function via the usual rela- 
tionship 

(n (r) n (V) ; t) = n8 (r - r') + n 2 g (r, r'; t) (9) 

where the brackets indicate an average over the (time-dependent) two-body distribution 
function. 
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B. Uniform Shear flow 



To induce shear flow, modified periodic boundary conditions are usedpq]. In the x- 
and z-directions, periodic boundaries are applied whereas in the y- direction, an atom with 
coordinates = (xj,yj,Zj) and momentum pj = (p X i,Pyi,Pzi) will have images with q = 
(xi + aLt,yi + L,Zi) and p = (p X i + aL,p yi ,p zi ) where t is the time and the parameter a, 
having the units of inverse time, is the shear rate. These are just periodic boundaries applied 
to the coordinates q£ = q« — atyix and p^ = p» — ayix which are the atomic coordinates in 
the local rest frame of a system undergoing uniform shear flow and are the standard means 
by which such a flow is simulated. 



It is known |]10| that this combination of dynamics and boundary conditions allows for an 
exact solution of the macroscopic conservation laws in which the local density is constant, 
the local flow velocity is given by v(r) = ar y x and the local temperature, defined as the 
excess kinetic energy relative to the flow field, is spatially uniform in the co-moving frame 
and increase as 

-nk B -T = -aP xy + F (10) 

where is the macroscopic pressure tensor, which is also spatially uniform, and the last 
term on the right represents the effect of the external forces. Typically, an external force, or 
thermostat, is included such that the right hand side of this equation vanishes, thus giving 
a constant temperature and allowing for the possibility of a stationary state. Here, it is 
assumed that for shear rates below the ordering transition, all thermostats are equivalent |2B 
27]. 



The one-body distribution function of a sheared and thermostated fluid of hard-spheres 
has been studied in considerable detail|l(|,[]ll|,and may be approximated as 



/(q,p) = P[^) [det(A)]- 1/2 exp(--/3^;A i , 1 ) (11) 

= 5 i:i + A i3 (12) 

where the (constant) matrix of coefficients is defined implicitly as the solution of 

Q> ($xi$jy + fiyifijx + ^xi-Ajy + 5 x jAiy) = C^j + ] m -Aj m (13) 
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with 



=f*J dq5( q -l)q r q J ^w 2 e- w2 / 4 + ±(w 2 + 2)4,(w)} 



(14) 



c ~ w2/4 + 6iP(w) \ q mi q m + { ^=e- w2 / 4 + i/;(w) \ A ijlm 



7T I V/7T 



CgL = ^PX y rfq - 1) 

+ S im qjqi + Sjiqiqmidjrnqiqi) 
ip{w) = w ^erf — lj 
w = aq x q y 

together with the condition that Tr(A) = 0. This approximation gives a semi-quantitative 
description of effects of strong shear such as shear thinning and normal stresses as well as 
being positive-definite at all shear rates. The pdf at contact within this approximation is 
found to be 

5 (r 12 - a) g{ Tl , r 2 ) = 5 (r 12 - a) Xo (l - erf (-L = ) ) (l 5 ) 

V V 2cx V 1 + Ai m r m r l2m J J 

which follows directly from Eq. (|l3|) and (P). Here and below, Xo is taken to be the equilibrium 
value of the pdf at contact as calculated in the Carnahan-Starling approximation. From this, 
the projections of the pdf at contact onto the spherical harmonics may be calculated as 

M lm = J dr 12 Y* m (r 12 ) 5 (r 12 - a) g{r x , r 2 ). (16) 

Note that from this point, the dependence of all quantities on time is being suppressed since 
we work in a steady state. 

III. THE NONEQUILIBRIUM PAIR DISTRIBUTION FUNCTION 

The Enskog approximation gives information about quantities at contact. In order to un- 
derstand the pdf for finite separations, several different approaches have been suggested. Two 
well-known proposals, the kinetic model of Hess[^] and the fluctuation model of Ronis|T3 



involve phenomenological parameters which can be fixed by requiring that they reproduce 
one of the moments M\ m as calculated from [IB]. These theories share the property that in 
equilibrium, they reduce to the equilibrium pdf so that it is reasonable to attempt to repro- 
duce local information such as the moments at contact. This contrasts with calculations of 



the density autocorrelation function based on kinetic theory |Tj|,@ or, equivalently, long- 



wavelength Langevin models[18] which can only give information at asymptotically large 



separations for which the information at contact is not relevant. 
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A. Kinetic Model 



The second equation of the BBGKY hierarchy is 



d ( d d d \ — 

—p 2 (x 1 ,x 2 ;t) + 22 \ Pi . ' ^7 + <li " *o* • 7^ + 7^7 • F(xi)J p 2 (x 1 ,x 2 ;t) + T-(12)p 2 (xi,x 2 ;t) 

-N J d3 (T_(13) + T_(23)) p^x^x^t) (17) 



and integrating over the momenta gives an equation for the pdf 

d . . , . d . . 



dpidp2T_(12)p 2 (xi,X2;t) -n J d Pl dp 2 J dS (T_(13) + T_(23)) p 3 (xi, x 2 ; $18) 



The kinetic models studied by Hess et al||12|| consist in replacing the complicated right hand 
side of this equation my a simpler diffusion or relaxation model constrained only by the 
requirement that it force a relaxation towards the equilibrium state. The simplest model 
then takes the form 

d d 

7^#2(qi2; t) + qi2 ■ *a* ■ -^—g 2 (sb2\ t) = -T (g 2 (q 12 ; t) - g 2 q (qi 2 )) (19) 
or, Fourier transforming, 

-^ 2 (k 12 ; t) - k 12 • *a> ■ g^U^m t) = -r (h(k 12] t) - h e 2 %k 12 )) . (20) 
The solution to Eq.([BJ), under the assumption of stationarity, is 

poo 

<? 2 (q 12 ) = / rf 7 e-^(g 12 (a 7 r)) (21) 
Jo 

where 

qi2(oT) = (Qi2x ~ aiqi2y, quy, q\2z) ■ (22) 
In Fourier space, this becomes the solution to Eq.(^C|) 



/>oo 

Jo 



with 

ki 2 (-a7r) = (ki2x, k\2y + a^Tk 12x , k 12z ) . (23) 

As alluded to above, the relaxation time appearing in this model can be fixed by requiring 
that the model reproduce one of the moments Mi m . 



B. Langevin Model for Density Fluctuations 



The Langevin model of Ronis[|13[] consists of a convective-diffusion equation for the decay 
of density fluctuations which, in Fourier space, appears as 

JUn (k, t) - k ■ <a* T ■ (k, t) — D (k) k 2 5n (k, t) = f (k, t) (24) 

where D (k) is a wave- vector-dependent diffusion constant and / (k, t) is a fluctuating force 
representing the neglected degrees of freedom. The fluctuating force is approximated as 
delta- function correlated in both wave- vector and time, with amplitude D{k)So{k)k 2 where 
So(k) is the equilibrium static structure factor. The pdf is obtained by solving for the 
density fluctuation as a functional of the force and evaluation of the equal-time density- 
density correlation function with the result 

~ r°° / n 

h(k 12 ) = / <*7 D(k(-a-f))h (k(-a-f))k 2 (-a-f) exp I- d-f'D (Jfe(-aY)) k 2 {-ai) 

(25) 

so that it is seen that the particular choice for the autocorrelation of the force leads to the 
correct result in equilibrium. The similarity between this and the Hess' model is apparent 
and in fact the same result is obtained if the relaxation time in the latter, Eq.(^), is 
taken to be wave-vector-dependent with T(k) = D(k)k 2 . To close the model, Ronis uses 
D{k) = D Q /So(k) with D a constant which he takes to be the equilibrium self-diffusion 
constant although in the present circumstances, it will be fixed by the requirement that the 
model give the correct value of M22. Finally, the nonequilibrium correction can be written 
more explicitly by means of an integration by parts which gives 

h{k l2 ) -h {k 12 ) = a Hdj k -4r-^p- K{K-^l)) exp ( - P diD (k(-aj')) k 2 {-ai\ 



H-ai) uv v UJ r V Jo 



(26) 



where h' {k) = ^h {k). The same result has recently been derived PS[ using a random phase 
approximation in the context of a Langevin model for the atomic coordinates. 
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C. Nonequilibrium GMSA 



As in equilibrium, define the direct correlation function as usual through the Ornstein- 
Zernike equation 

h(r u r 2 ) = c(r h r 2 ) + J dr 3 c(ri, r 3 )p(r 3 )/i(r 3 , r 2 ) (27) 

where h(ri, r 2 ) = g{ri, r 2 ) — 1 is the structure function. The pdf must satisfy the boundary 
condition that the probability for two atoms to interpenetrate is zero so that 

/i(r a ,r 2 )=-l | ri -r 2 |<a (28) 

The OZ equation can be solved for both the structure function and the direct correlation 
function provided that this is supplemented by a closure relation between the two. 

In equilibrium, the relation between the pdf and the direct correlation function may be 
written as 

c(r h r 2 ) = higfa, r 2 ) - h(r 1 , r 2 ) + u(n, r 2 ) + B(r 1 , r 2 ) (29) 

where v(ri, r 2 ) is the pair potential and B(r\, r 2 ) is the bridge function which is not generally 
known in closed form. The integral equation approach to liquid state structure can be written 
in terms of various approximations to the bridge function. Setting B = yields the hyper- 
netted chain equation and further approximating In g = In [1 + h] ~ h, or B = h — In (g), 
yields the Percus-Yevik approximation. A number of other approximations exist, including 
schemes such as that of Rogers and Young |29|| and of Zerah and Hansen 0], which involve 
a parameterization of the right hand side of Eq.(^). For hard-core potentials, the Percus- 
Yevik approximation reduces to the statement c(ri, r 2 ) = for |ri — r 2 | > o and the GMSA 
replaces the right hand side by a Yukawa function with parameters adjusted to give a known 
equation of state. In this case, these have been shown to be the first two steps in a systematic 



expansion of the tail of the dcf with little underlying physical approximation f20|. [2T| . [2^ , |23 
It thus becomes natural to carry over this model to the nonequilibrium state so that the 
closure condition is expressed in terms of a similar parameterization of the tail of the dcf 
giving the form 

c(ri,r 2 ) = } j A i K i (r 1 ,r 2 ) |r 1 -r 2 |>a (30) 
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for some set of basis functions {Ki(ri, r 2 )}. As it stands, Eq.(pO[) is quite general and the 
physical approximation will be to truncate and parameterize this expansion as discussed 
below. 

The problem of solving the OZ equation is now formally equivalent to that of the case of 
molecular fluids and similar techniques can be used|JT|. To begin, one expands the angular 
dependence of the dcf, the pdf and the boundary conditions in terms of spherical harmonics 
so that 

hfa, r 2 ; t) = ^ hi m (r 12 ; t)Y lm (r 12 ) (31) 

lm 

with similar expansions for the other quantities. Using Rayleigh's expansion of a plane wave 
in terms of radial and angular functions and the addition theorem for spherical harmonics, 
it is easy to show [[31]] that the Fourier transform of such an expansion has the form 

h(k u k 2 ;t) = (2n) 3 5(k 1 + k 2 )J2hi m (h;t)Y lm (ki) (32) 

lm 

with the coefficients defined in terms of Hankel transforms 



h lm (k; t) = Am 1 / r 2 dr ji(kr)h lm (r) . (33) 
Jo 

The Fourier transform of the OZ equation then becomes 

hi m (k) = Ci m (k)+nc Vm ,(k)hi>' m ''(k) J dk ^/ m /(k)^// m »(-k)Y^(k) (34) 

1 V 
(k) + n—= ^ ^2 j4 (M'i'>,m')QwW'iiw(*:) 



'An 

' \l'-l"\<l<l'+l" m'=-V 



where 



A(l, r, l", m, to') = (-l) v+2m J (2// + 1)(2 ^ + 1) C(r, l", l\0Q0)C(l', I", l\m', m - m', to) 

V 2/ ~\- 1 

(35) 

and the last line is a well-known result which follows from the Wigner-Eckart theorem |32] . 



Equation fl3~4|) together with Eq.(^) and the exact condition, Eq.(^) serve to define the 
integral equation. 

In the theory of molecular liquids^], auxiliary functions are usually defined as 

f (r) = J ^^^dk^UUk) ^ven 
'"" ^Io°°k 2 dk^7 lm (k) /odd 
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where fi m (k) could be either hi m (k) or ci m (k). For even values of I, the only ones of con- 
cern below, this means that the Hankel-transform of the original functions is the Fourier- 
transform of the auxiliary function. This is useful for numerical work but more important 
is that the auxiliary functions tend to be of shorter range than the original functions. To 
see this, we need the relation between the auxiliary functions and original functions in real 
space 

r°° 1 

fL(r) = fUr) - / r> 2 dr> —P[ {r/r>) f lm {r>) (37) 

J r 

where P((u) is the derivative of the I— th Legendre polynomial with respect to its argument. 
Direct calculation shows that if fi m (r) = r~ n then // m (r) = provided that 2 < n < I + 3 
. This transformation therefore removes a subset of long-ranged decays and is consequently 
degenerate. The inverse transformation is 

[¥] 

fUr) = (~l) l fL(r) - (-1)' / r'Hr' ^P[ (/ '/r) f lm {r>) + £ A lm , n r-^ +1 \ (38) 

JU n=l 

where the relationship holds for any choice of the coefficients Ai m>n - an expression of the 
degeneracy of the original transformation. In the present application, since h\ m (0) is, from 
eq . fl28D , we will always have Ai m ^ n = whereas for the dcf, no such statement can be made. 

Further discussion of the details of the solution of these equations is given in appendix 
^ and only some of the conclusions of that analysis are stated here. First, because of 
the symmetry of the boundary conditions, only coefficients corresponding to even values 
of I and m will be nonzero. Second, one expects that, because of the symmetry of the 
flow, the dominant nonequilibrium contributions come from I = 2 and m = ±2 (since 
Y 2 2 + Y 2 -2 xy) and indeed this can be verified for the pdf at contact by calculations using 
Eq.flI5p. In the case that we keep only these contributions to the OZ equation, as well as the 
I — m = component necessary to describe the equilibrium contribution, the problem can be 
transformed into the solution of two one-dimensional OZ equations with the Yukawa closure 
and an analytic solution is possible as described in appendix |A|. This approximation should 
be understood in the spirit of a truncation of a moment solution rather than an expansion 
in the shear rate and in fact, all of the expressions presented below depend nonlinearly on 
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the shear rate. The result is that 



h(r) = /i o(r)r o(r)+2/i 22ir (r)Rer 22 (r)-2/ i22ii (r)Imr 22 (r) 
where the angular average of the pdf is given by 



(39) 



with n± 



n 



1 ± Id 



h 00 (r) 

hp* ( B - 



In 



(n + h + {r\ n + ) + n_/i_(r; nJ)) 



(40) 



22,1/ 



and where the constant, -B 22 ,i, is related to A 22 ,\ oc- 
curing in eqfl3"8|). As discussed in detail in the appendix, the functions h±(r;n) may be 
expressed in terms of the solution of the OZ equation for homogeneous system (e.g., the 
Percus-Yevik solution if the right hand side of Eq.(|30D is set to zero or the known analytic 
solution of the N- Yukawa closure^] if the basis functions are Yukawas) for density n±. The 
anisotropic component is given by 

£?22,1 



lm.(h' 22 (r;B 22 ,i)) 



6 (r - 1) 



h' 22 (r';t) + 



4 r r' 2 dr' ti 



r 3 ./ 1 



22 



(r': B 



22,l y 



(41) 



(n + h + (r; n + , ) — n_/i_(r; n_)) 



2 \a\ n 

Re(h' 22 (r)) = Im(/ i , 22 (r))Re(M 22 )/Im(M 22 ) 

Aside from the truncation of the OZ hierarchy, these results are independent of any assump- 
tion about the closure condition given in Eq. (|30|) . 

In equilibrium, the GMSA is based on the choice of a Yukawa function for the tail of the 
dcf. This is motivated by the expectation that the tail is short-ranged and, then, because 
a Yukawa closure is analytically tractable. As discussed above, recent work has shown that 
the Yukawa may be thought of as the first term in a systematic expansion thus removing 
some of the arbitrarity of this choice. In the same spirit, we therefore take as the principle 
hypothesis of the extension of the GMSA to nonequilibrium systems that the tail of the 
auxiliary dcf function can be expanded as 



6 (r i2 - a) c'(ri, r 2 ) = ^ v' lm (r 12 )Yi m (r 12 ) 



with 



lm 



Lm 



r) ~ -Ki m exp(-zi m r) 
r 



(42) 



(43) 
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For the spherically symmetric component I = m = this is just the Yukawa closure as in 
equilibrium. Because of the non-uniqueness of the relation between the dcf and the auxiliary 
def, see the discussion after Eg. (p8[) , this corresponds to a closure of the actual dcf of the 
form 



6 (r 12 - cr) c, TO (n, r 2 ) = v lm (r 12 ) + ^ A lm ^ 2n+1 \ (44) 

n=l 

Since the dcf is of no significance in the present context, the values of Ai m ^ n are left indeter- 
minate at this stage. 

As formulated, the truncated analytic model has 5 parameters corresponding to the 
amplitudes and length scales of the two Yukawa terms and the constant -B 22j i appearing in 
Eqs.(]39|)-pTD. These parameters are constrained by two boundary conditions consisting of 
the values of M 0Q and M 22 . In the first application to shear flow||, the model was simplified 
by setting Kqq = K 22 = 0. The justification for this was simplicity, since there is then 
only the non- uniqueness parameter, -B 22j i, and it was shown that a value could be found 
which simultaneously satisfied both boundary conditions reasonably well. However, recent 
estimates @] indicate that h 22 (r) decays faster than 1/r 3 for large r leading to the condition 
that the inverse-cube terms in Eq. (|4~l~D vanish in the limit r — > 00 giving the constraint 



00 

S 22jl = 3 J r' 2 dr' h' 22 (r';B 22jl ). (45) 

This still leaves two parameters undetermined. One of these will be eliminated by taking 
the length scale of the length scale of the v' 00 (r) function to be fixed at its equilibrium value 
and to only allow the amplitude to be adjusted so as to reproduce M 0Q . This still leaves one 
undetermined parameter which can be taken to be z 22 . In an application of this approach 
to granular fluids, a similar indeterminacy was solved by insisting that the compressibility 
equation continue to hold in the nonequilibrium state. Here, this is not useful because v 22 (r) 
has little influence on /ioo which would be the object fixed by such a relation. (In fact, this 
could be used as an alternate means of fixing z Q0 .) With no other exact or well-approximated 
property to fit, it seems appropriate to try to minimize the perturbation of the tail of the 
dcf. The tail of the full dcf is found to be 

/ \ v exp (-222 (r - 1)) r 2 z\ 2 - 3z 22 r - 3 3 f°° 2 

c 22 (r) = K 22 5-= / c 22 (x)x z dx for r > a (46) 

r r z z 22 H Jo 
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which is clearly short ranged if the last term on the right vanishes as it in fact does in the 
present approximation as a result of the condition given in Eq.(^). For large separations, 
the tail is therefore the same as that of the auxiliary function, a Yukawa and this cannot be 
changed by any condition on z 2 2- At short range, the Yukawa is modified and one possibility 
that suggests itself is to demand that at contact, the tail assume its equilibrium value - 
namely, zero. This implies that Z22 = § + Iv^T ~ 3. 8 which is the value used below. 



IV. MOLECULAR DYNAMICS SIMULATIONS 



In order to evaluate the model for the structure proposed above, I have performed molecu- 
lar dynamics simulations of sheared hard spheres in three dimensions. The simulation makes 
use of Lees-Edwards boundary conditions to impose the shear and the heating is controlled 
by periodically rescaling the velocities. Specifically, to maintain an average temperature 
To, the velocities are rescaled to give an instantaneous temperature of 0.95Xo whenever the 
instantaneous temperature exceeds 1.05To. The simulations reported here were performed 
using 500 atoms except where noted below. In all cubic simulation cell was used. 

The equilibration procedure consisted of first creating an equilibrium liquid at the desired 
density. After 10 7 collisions, the shear rate was then instantaneously set to the desired value 
and the system allowed to relax under the Lees-Edwards boundary conditions for another 
10 7 collisions. Finally, the simulation was extended for another 10 7 collisions during which 
statistical averages were accumulated under the ergodic hypothesis. In order to estimate 



the accuracy of the quantities obtained, Erpenbeck's pooling method |35[ was used whereby 
averages were accumulated over periods of 10 5 collisions and stored. The reported values 
were subsequently computed by averaging these partial-averages and the standard error of 
these partial averages, i.e. the standard deviation divided by the square root of the number 
of observations, used as an estimate of their accuracy. Except as noted below, the error 
bars in all figures are smaller than the symbols used to display the data. Finally, all simu- 
lations were performed for reduced densities of n* = na 3 = 0.1, 0.25, 0.5 and 0.75. Based on 
the difference of the equilibrium pdf at contact from its low density value, namely Xo — 1> 
these densities correspond to low density (xo — 1-14), moderately dense (xo = 1-4), dense 
(Xo = 2.15) and very dense (xo = 3) fluids respectively. Shear rates are reported in units of 
the Boltzmann collision time a* = a/ (^4n* y/TrksT) . 



16 



20 



15 



S 10 



04 



•- 



0.5 



1.5 



FIG. 1: i?e(Moo) as a function of the reduced shear rate, a* for densities of 0.1 (circles), 0.25 
(squares), 0.5 (diamonds) and 0.75 (triangles). The lines are the predicted values calculated from 
eq.(|i~5|),Note the non-monotonic behavior at the highest density. 

A. The pair distribution at contact 

Figures |I],|2] and show the projections M; m of the pdf at contact onto the spherical 
harmonics for Im = 00, 22 and 44 as a function of shear rate which accounts for angular 
dependencies of the form 1, q x q y and q^—q^ respectively. The spherical average of the pdf, goo, 
shows little variation with shear rate for n* = 0.1 and 0.25 and only begins to show significant 
variation above a* = 0.5 for n* = 0.5 while at the highest density, significant variation is 
observed for all shear rates and, unlike at lower density, the curve is not monotonic. In all 
cases, the Enskog prediction is a slight decrease with increasing shear rate which is confirmed 
in the low density data. The system at intermediate density is consistent with the model for 
small shear rates but shows an increase at higher shear rates as does the high density system 
at all shear rates, in qualitative disagreement with the model. The major nonequilibrium 
contribution to the structure resides in the Im (522) components which show qualitatively 
similar behavior: agreement with the model at low density and all shear rates and for 
low shear rates at moderate density with significant disagreement at moderate density and 
high shear rates and at all shear rates at high density. The next largest nonequilibrium 
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FIG. 2: Same as Fig. [T] but showing M 2 2- 
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FIG. 3: Same as Fig. |l| but showing M44. 

contribution, Re (#44), is seen to be nearly an order of magnitude smaller than Im (#22) and 
it is also poorly described by the model. 

These results show that the largest nonequilibrium contributions to the pdf occur in the 
four directions (±75, ±7?) 0). To give a direct overview of the accuracy of the models, the 
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FIG. 4: The value of \ (gg^ + 910,2) as defined in eq.47. The labeling is the same as in Fig. |T|. 
average of the pdf at contact over a number of angular bins, defined as 

/•- l+(m+l)6 m /-(n+l)^ 

9mn = d cos 6 I d(f) g (q) , (47) 

J — 1+mSx J «<5</> 

was monitored during the simulations for 5 X = i and 5^ = f^. Figure |] shows a comparison 
between theory and simulation of | (#9,2 + #10,2), e.g. of the pair distribution at contact 
averaged over the area —0.1 < x < 0.1 and 0.2tt < <fi < 0.37T, for the various densities. This 
patch is centered on the direction (^, ^,0) for which the deviations from equilibrium are 
largest. The most striking feature of these results is that the pdf drops with increasing shear 
until it becomes vanishingly small indicating that no collisions take place in that direction. 
The figure also shows the predicted values based on the generalized assumption of molecular 
chaos, Eq.flTEl), averaged (numerically) over the same solid angle. It is evident that the 
model works quite well at the lowest densities, is reasonable at the intermediate density and 
is only qualitatively correct at the highest density. In order to visualize the full directional 
variation of the pdf at contact, Figure |5] shows the spatial variation of the pdf averaged over 
the same sized solid angle for the whole range of values of from — 71 to 7r for fixed shear 
rates of 1.0 for n* = 0.1,0.25 and Fig. ^| shows the same for a* = 0.6 and n* = 0.5,0.75. 
(The reason for choosing a lower shear for the higher densities will become apparent below.) 
For all but the highest densities, the spatial variation is consistent with the model, Eq. (|15l), 
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FIG. 5: The pdf at contact averaged over a solid angle denned by —0.1 < x < 0.1 and ^ < ^ < 



m+l 
10 



for —10 < m < 10. The shear rate is fixed at a* = 1.0 and the information for n* = 0.1, 0.25 



is labeled as in Fig.|]. 

whereas for n* = 0.75 the agreement is poor. Indeed, the simulation data in the latter case 
is erratic and appears to be a superposition of a periodic function with spikes near = 
and = 7r which corresponds to the directions parallel and antiparallel to the flow (i.e., ±x). 
Figure |7| shows the same quantity for n* = 0.5 and a* = 1.0 and the same superposition of 
features is apparent. There are three possible causes for deviations from the model : shear- 
induced ordering, inaccuracy of the one-body distribution and breakdown of the assumption 
of molecular chaos. The structural anomalies at high shear rate and high density suggest 
the former. 



It has been known for some time that hard spheres undergo an ordering transition at high 
shear rates|3C|. The nature of the ordered phase remains uncertain and appears to depend on 
the type of thermostat used^7|. For the simple rescaling thermostat used here, an ordering 
first into plane perpendicular to the direction of the gradient (here, the y-direction) and 
then into strings oriented along the direction of flow, and in a hexagonal pattern in the 
plane perpendicular to the flow, has been reported [35]. As a quantitative measure of such 
an ordering, the average density in a tube, oriented along the direction of flow has been 
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FIG. 6: The same as Fig. | for a* = 0.6 and n* = 0.5, 0.75. 
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FIG. 7: The same as Fig. | for a* = 1.0 and n* = 0.5. 



monitored 101. This is defined as 



t{u) 



niru 2 L \ N 



21 




FIG. 8: The tube density t(a) as a function of the shear rate. The labeling is the same as in Fig. 
|] except that here the lines are only a guide to the eye with the baseline, t(a) = indicated by the 
thick line. Open symbols are from simulations of 108 atoms. 

which can be written in terms of the pdf as 

t{u) = -L- fdt g(r)Q {u 2 -y 2 - z 2 ) (49) 

7TU L J 

= ^ + ^ + ^1^ %)6 (r 2 -1)0 (u 2 -y 2 - z 2 ) 

and gives the average density, relative to the bulk, observed along a tube of length L and 
radius u centered on an atom. In the limit of large L, the last term on the right will only give 
a non-zero contribution if long-range correlations in the direction of the flow are present (as 
they would be for a " string phase" ) so that any deviation from the equilibrium value could be 
attributed to the formation of such correlations. However, for the small systems considered 
here, the last term will give a nonzero contribution in all circumstances and variations of 
the tube density with the shear rate could be due to variations in the pdf which nevertheless 
do not involve long-ranged correlations. Figure |8] therefore shows the tube density as a 
function of shear rate for systems of both 108 atoms and 500 atoms (giving L = 6 and 10 
respectively). For the lowest densities, the tube density actually decreases with increasing 
shear rate with the decrease being larger for the smaller system. Noting that the size of 
the effect is roughly in inverse proportion to the length of the systems and independent of 
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the density, this would appear to be primarily a finite size effect leading to the conclusion 
that in neither case is there evidence of shear-induced changes in the density in the infinite 
system limit. For n* = 0.5, the tube density is roughly constant until above a* = 0.6 at 
which point it begins to rise steadily for both system sizes. Although the magnitude of the 
increase is also a function of system size, the relative variation between the two systems is 
much less than the over-all increase leading to the conclusion that the increase is due to 
the development of long-range order. This is consistent with our previously reported results 
indicating that shear- induced ordering takes place at this point 0. Finally, at high density, 
the tube density increases dramatically above a* = 0.2 indicating an ordering transition 
at that point. This behavior, in this case nearly independent of system size, supports the 
conclusion that the deviation of the data from the molecular chaos hypothesis is due to 
shear- induced ordering for n* = 0.5 and a* > 0.6 and for virtually all of the data for the 
high density system. It is also consistent with the structural data which show spikes in the 
pdf corresponding to an increase in collisions in the direction of the flow. 

We are then only left with the poor agreement of the model for the #44 to explain. It 
seems likely that this is simply due to the inadequacy of the information supplied for the 
one-body distribution. Since the distribution is accurate only up to second moments of the 
velocity, it is reasonable that the calculation is only accurate up to second order in the 
unit vectors. Since these contributions are in any case small compared to the dominant #22 
terms in the region of validity of the model, this aspect of the problem has not been pursued 
further. 



B. The pair distribution function at finite separations 

Here, attention is restricted to the domain of densities and shear rates below the ordering 
transition. In general, the components gi m (r) were estimated during the simulations by 
evaluating 

9i m (r) = -^— T^E (^- r ) ( r + rfr -^)^(^) ( 5 °) 

iv samples , 1 111 ... 

c samples i<] 

where the inner sum is an estimate of gi m (r) based on a snapshot of the system and the 
outer sum indicates an average over many different snapshots. The results presented here are 
based on snapshots taken every 10,000 collisions. Ideally, one would like to replace the outer 
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FIG. 9: The functions Re (goo(r)) /v47r, upper curves, and Im (g22( r )), lower curves, for n* = 0.1 
and a* = 1.0 from simulation(circles), the GMSA with the Yukawa closure (full lines) and the 
GMSA with the tail of the dcf set to zero (dotted lines). 

sum by a continuous time-average in the spirit of the ergodic theorem but the computational 
expense would be prohibitive. 

Figures |9l,|T0| and |TT] show Re(g 00 (r)) and -/m(g 2 2( r )) for n * = 0.1,0.25 and 0.5 and 
a* = 1.0, 1.0 and 0.6 respectively as determined from the simulations. While the spher- 
ically symmetric component is little changed from equilibrium, the main nonequilibrium 
component, IiTL(g 2 2{r)), is comparable in magnitude near the core to the equilibrium pdf 
but decays rapidly and is difficult to measure beyond about two hard sphere diameters. Also 
shown in these figures are the results of the GMSA model with both the Yukawa closure 
described above and the simpler closure in which the full dcf is set equal to zero outside 
the core. Both give a good description of the main features of the structure including the 
location of the sign changes and the amplitude and wavelength of the oscillations in g 2 2 
with the Yukawa closure being obviously superior in all cases. It is interesting that the 
largest discrepancy occur for the lowest density Based on the observed fluctuations in the 
identical determination of gi m {f) for this density in equilibrium, the statistical errors for this 
system appear to account for no more that half the deviation seen in Fig.|| One noticeable 
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FIG. 10: Same as Fig.B for n* = 0.25 and a* = 1.0. 




FIG. 11: Same as Fig. Q for n* = 0.5 and a* 



0.6. 



characteristic of the simulation data for n* = 0.1 is that Im(g 2 2) appears to be systemat- 
ically below the GMSA model and indeed systematically below zero away from the core. 
This suggests that the deviation might be a sign of the slow, algebraic decay predicted by 



mode-coupling theory [34] and which is not incorporated in the GMSA model. 
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FIG. 12: The functions Re(goo{r)) j\f^n and Im (<?22(?")) for n* = 0.5 and a* = 0.6 from simula- 
tion(squares and circles, respectively), the Hess model with the small er parameter (full and dotted 
lines, respectively) and with the larger parameter (dashed and dot-dash lines, respectively). 

Figures ^ and [13] show the same simulation data as Fig.|T^ together with the numerical 
evaluation of the models of Hess (eq. (|B^)) and Ronis (eq. (^§)) performed by means of a Monte 
Carlo integration using the VEGAS algorithm ^ , [38j j39[ . Because it can be formulated in 
real space, the Hess model requires one fewer integral than the Ronis model so that the 
evaluations are quicker and more accurate. In both cases, the parameters were adjusted so 
as to reproduce the calculated value of 7m(M 2 2). Curiously, there are two values of the 
free parameter in the Hess model that satisfy this constraint. For the smaller of the two 
values, Re(goo(r)) is almost unchanged from equilibrium and Im(g22{i r )) is only nonzero in a 
very small region near the core. Both components are non-zero in a small region inside the 
core. The larger value of the parameter gives rise to a substantial deviation in i?e(goo( r )) 
which is therefore not in good agreement with the simulation data. In contrast, Im(g22(r)) 
is in qualitative agreement with the data outside the core. In this case, both components 
take on substantial values inside the core. The Ronis model, shown in Fig. is similar 
to the large-parameter version of the Hess model. The spherically symmetric component 
is modeled somewhat better but Im{g22{ r )) is somewhat worse than with the Hess model. 
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FIG. 13: The functions Re(goo{r)) j\f^n and Im (<?22(?")) for n* = 0.5 and a* = 0.6 from simula- 
tion(squares and circles, respectively), and the Ronis model (full and dotted lines, respectively). 
The error bars are the standard errors reported by the VEGAS algorithm used to evaluate the 
model. 

The behavior inside the core is also worse, with Re(goo(r)) even taking on negative values. 



V. DISCUSSION 

The goal of the work presented here has been to describe the density-density correlation 
function in a sheared hard-sphere fluid over a range of densities, length scales and shear rates. 
It was shown that the Enskog approximation for the velocity correlations provides a basis 
for calculating the the density-density correlation function at contact and that the results 
hold up well when compared to simulation, even for conditions of moderate density and high 
shear rates. Indeed, deviations from the Enskog predictions are primarily attributable to 
the high-shear rate phase transition and, as noted previously!!, appear to signal its onset. 
For finite separations, the nonequilibrium GMSA was shown to provide a framework within 
which the known information about the correlation function at contact, coming from the 
Enskog approximation could be used to provide an accurate description of the dominant 
effects at finite separations. While power-law decays arise naturally in the solution of the 
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anisotropic OZ equation, the inverse-cubic decay in g 22 (r) is in contradiction to the l/r 17//3 
decay predicted by recent mode-coupling calculations [[34]]. Furthermore, those calculations 
indicate a 1/r 11 / 3 decay in goo( r ) f° r which there is no analog in the OZ solution. Of 
course, such algebraic decays could find their origin in the auxiliary functions and, indeed, 
might be used as boundary conditions in place of the Yukawa closure of the GMSA. On 
the other hand, the mode-coupling results are the result of a number of approximations 
(linearized Navier-Stokes-Langevin model solved perturbatively) and may not be giving the 
true asymptotic behavior. Altogether, the question of the actual form of these algebraic 
decays must be considered to be unresolved at this point since the data presented here is 
adequately fitted without them (except, possibly, for n* = 0.1). More extensive simulations 
which can provide better statistics for the decay at separations significantly above two or 
three hard-sphere diameters would be required in order to resolve this issue. What can be 
said here is that even if algebraic tails exist, they must be significantly weaker than the 
dramatic nonequilibrium contributions seen near the core. 

Two other, closely related, theories for the pdf were also considered and shown to capture 
some of the qualitative behavior of the pdf but both suffer from the unphysical prediction 
of nonzero probabilities inside the core. Before dismissing this class of theory on this basis, 
it is interesting to consider whether the main failing could be eliminated. In order to show 
that this is indeed the case, consider the second equation of the BBGKY hierarchy for hard 
spheres 



p 2 (xi,x 2 ;t) (51) 



l+£(p^ + *- v -^ + lr F H) +r - (12) _ 

= -N Jd3 (T_(13)+T_(23))p 3 (x 1 ,x 2 ,x a ;t) 

First, observe that in the most general case the distribution must have the form 

p 2 (x h x 2 ; t) = O (ri2 - a) p 2 (xt,x 2 ; t) (52) 

since the atoms cannot interpenetrate. Substituting this into Eq. flBTD , one finds two terms 
proportional to 5 (r 12 — a) : the first coming from the action of the spatial gradients on the 
step-function in eq.([5^) and the second from the collisional operator T_(12). These must 
cancel so that their coefficients must be equal and this just gives the relation between pre- 
and post-collisional distribution functions used to derive Eq.(|). The remaining regular part 
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of the kinetic equation then reads 
G (qi2 - cr 



i=l,2 v 



0(h <h\, dp't 

-NQ(q 12 -a) J d3 (T_(13) +T^23))p z (x 1 ,x 2 ,x 3 ;t) 



p 2 (xi,x 2 ;t) (53) 



where an analogous decomposition of the 3-body distribution has been introduced. Inte- 
grating over all momenta and discarding surface terms then yields 



6 (q 12 - cr 
-0 (q 12 - 



d 



d 



5qi2 



2/(qaz;*) 



(54) 



a) n 



J dpirfp 2 J d3 (T_(13) +T_(23)) p 3 (x 1 ,x 2 ,x 3 ;t) 



where y(qi,q2;t) is the nonequilibrium cavity function and is related to the pdf by 
g(qi,q2;t) = O (gi2 — c) y(qi, q 2 ; t). This suggests that a more physical approximation 
would be to make a relaxation approximation for the cavity function of the form 



6 
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9 , 

^y(qi2; t) + qi2 



d 



l/(qi2! t) - / dr A(q X2 - r) (y(r; t) - y (r)) 



= 0. 

(55) 



There is no reason at this point to keep the step-function in this equation since any solution 
valid for all separations will be valid outside the core. Then, this gives in Fourier space 



J^(k; t) - k ■ V T ■ J-y(k; t) = A(k) (y(k; t) - gb(k)) 



(56) 



and for steady-state shear flow the solution is 

y(k) = y (k) + Tdtf ak ^- at) %(k (-at)) ) exp 



dt' A(k(-at')) 



(57) 



An extension of the Hess model would be to take A (k) = ^2i m Ai m Yi m (k) for some set of 
constants Ai m adjusted to give the correct moments at contact. A similar extension of the 
Ronis model is also possible. An investigation of these models will be left to a future study. 

In conclusion, it has been shown that the Enskog approximation for the pair correlations 
at contact, together with the GMSA model, provides a good description of the density 
autocorrelation function in a sheared fluid. The same techniques also give a good description 
of the pdf in granular fluids (modeled as inelastic hard spheres) 0] giving evidence that the 
approach is applicable to a variety of nonequilibrium systems. In both cases, simply knowing 
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the value of the pdf at contact, from the Enskog approximation, and applying the standard 
formalism of liquid-state theory are enough to give a description of features of the system 
arising solely from the nonequilibrium state. Further work will include the extension of this 
model to the description of static correlation functions involving temperature and velocity 
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APPENDIX A: SOLVING THE OZ MODEL FOR USF 



The OZ equations written in terms of the auxiliary functions are 

1 J. 

Km (*) = Am (*) + n-= E A Q> ™> (*) h 'l"™-m> (*) (Al) 

71 \l'-l"\<l<l'+l" m'=-V 

and the boundary conditions, which follow directly from Eq. (|37|) , are 

e(a - r)h' lm (r) = -V^S l0 S m0 + (1 - 5 l0 ) £ B lm , n r 2n (A2) 

n=0 

9 ( r - ff )4( r ) = ^m( r ) 

where the constant coefficients are functionals of the structure function 

1 f°° / r \ 

B ^T 2n = — dr ' E ^ + 3 ) P 2n + 1 ("J ^ ^ *) ' ( A3 ) 

n=0 ra=0 

Note that Eq. (|A3|) is not a self-consistency condition: rather, it will automatically be satisfied 
for any solution of eqs.(Al) and (|A^) as can be verified using the relation 

h lm (r; t) = ti lm (r; t) - \ £ (4n + 3) f r'dr' P 2n+1 ^ m (r' ; f) . (A4) 

r n=0 ^° V r / 

Instead, the significance of the constants is revealed by considering the equivalent relation 
for the dcf 

c lm (r; t) = c' lm (r; t) - \ J] (4n + 3) / r'dr' P 2 „+i f- j cj ro (r'; t) (A5) 

r n=o 70 V r / 
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which implies that 

6 (r - 1) q m (r; t) = 9 (r - 1) v'{ m (r; f) - (A6) 

l-i 

( r - !) 3 E ( 4n + 3 ) / r ' dr ' P ^+i (-) W« {r'\ t) - v' lm (r'; t)) 



where 



v'L(r) = «L (r) - 4 E ( 4n + 3 ) f r ' rfr ' P ^ + i (-) <i> (r) ■ (A7) 

r n=0 ^° V r / 

The reason for defining v" m (r) is precisely due to the non- uniqueness of this transformation, 

as discussed below Eq . (p8|) of the text. The point is that if we assumed a closure of the full 

f— 1 

dcf of the form Uj m (r) = (f>i m (r) + J2n=i ^im,n r ~^ 2n+1 \ then only 0/ m (r) would contribute 
to v' lm (r) and we would find that v" m (r) = 0/m( r )- We would then complete the problem by 
adjusting the constants B\ mn so that the boundary condition is satisfied, meaning 



E ^^" (2n+1) = 4 E ( 4n + 3 ) f r ' dr ' p < ' r 

n=l ra=0 ^° 



However, in the present application, we are not concerned about the full dcf and so the -B/ m , n 
are simply treated as free parameters. 

If we retain only the I = 0,m = and I = 2,m = ±2, components, the model can be 
reduced to the solution of a one dimensional OZ. The explicit form of the OZ equations in 
this case are 

Ko = c 00 + n \J^ c 00 h' oo + (c 22 h! 2 _ 2 + c 2 _ 2 h' 22 ^j (A8) 

^22 = C 22 + n ^4^ ( C 00^22 + C 22^0C>) 
h'2-2 = C2-2 + n j J~ (^00^2-2 + ^2-2^00 



with the boundary conditions 



Q(a - r)h' 00 (r) = 


-V47T 


Q(a - r)h' 2±2 (r) = 


— -02±2,O 


0(r - a)c' 00 (r) = 




9(r - o-)c' 2±2 (r) = 




limCta + e) = 


Moo 


lhn/i 2±2 (o- + e) = 


M 2±2 + B 
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(A9) 



First, note that hi^ m = {—l) m h* lm1 etc. which follows from the equivalent property of the 
spherical harmonics. Then, it is useful to separate the equations into their real and imaginary 
parts to get 



1 




Cqo^'oo + 2 



( C 22,r^22,r + C 22,J l 22,i^ 



(A10) 




^22,r _ C 22,r + n \j ^ (^ C 00^22,r + C 22,r^Oo) 



h 22,i = c 22 i + n\j — yc 00 li 22i + c 22 i n 00 j 



where h' 22 r = i?e ^/i' 22 j , etc. Second, because of the linearity in the 22-components of the 
last two equations and the boundary conditions on the core-values of the components of the 
structure function are constants, these equations are solved by taking 



^22,r — Xh' 2 2,i 



x = M 22)T /M 22)l 



(All) 



giving 



K 



00 — HlO 



Cnn+U 




? 00 ti 00 + 2(l + x 2 )? 22t hl 



22A 



1 



22,r — ^22, r 




c 99 r + n\l ^ [c 00 h 22 r + c 22 r h 00 



(A12) 



provided Re (v 22 ) = xlm (v 22 ) which we are free to impose. Now, define h(r; u) = h' 00 (r) + 
uh' 22i (r), which satisfies 

h(k, u) = c(k, u) + n \J~-^_ %A'oo + u {cQ Q h' 22 r + c 22 r h' 00 ) +2(1+ x 2 ) ^h' 22 J(A13) 
= c(k,u) + n\ -j—c(k,u)h(k,u) 

V 47T 



where the last line follows if and only if u — ±a/2 (1 + x 2 ). The boundary conditions are 
then 



©((7 - r)h(r; u) = —V^ — ulm(B 2 2,i) 
e(r-a)c(r,u) = (r) + uv 2 { 2 ] (r) 



(A14) 



Finally, introduce scaled quantities defined by h(r;u) = [y/An ± \u\ Im (-B 22 ,i)] H±(r) and 



32 



c(r;u) = [\/47r + \u\ Im (.822,1)] C±(r) which give 
H±(k) = C ± (k)+n±C±(k)H ± (k) 



(A15) 



S(a-r)H±(r) 
Q{r-a)C±{r) 



47r ± \u\ Im (B 2 2,i] 



Voo\r) ±\u\ 



An ± \u\ Im (-822,1) 



-1 -1 



^2 ( 2 0) ( r ) 



which can be recognized as the OZ equation for a (possibly negative) density n± 



n 



1 ± \u\ \ / j-ina. (-822,1) and some particular closure condition (so that this resembles 



the usual GMSA). For the simplest case, in which one takes v x 



00 1 



0, we have 



H+(r) 



' (0) 

hpy(r;n±) and the solution is trivial. If the tail functions,t> 00 (r) and v 



22 



are Yukawas, then we can use of the solution of Hoye and Blum for a closure consisting of 
a sum of Yukawas ||33| I . For completeness, we collect together the various transformations to 
see that 

VAtt 



Ko(r) 



2n 



2 \a\ n 
xh' 22 Ar) 



n + H + {r) + n_H_(r)) 
(n + H + {r) — n^H_(r)) 



(A16) 



. The value of Im (-822,1) is, of course, fixed by requiring that 

lim ti 22i (a + e) = -Im (-832,1) + Im (M 22 ) 

while the full pdf is 



h(r) = /4(r)y 00 (r) + 2h' 2 ' 2:r (r)R e Y 22 (r) - 2h / 2 \ l (r)lmY 22 (r) 



(A17) 



(A18) 



where 



^22 (0 



r' 2 dr' U 2 (r', r)h! 22 (r'; t) 

efr-i 



(A19) 



r' 2 dr' h! 22 (r'; t) 
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